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SUBCRITICAL l 1 '!. UTTER TESTING AND SYSTEM IDENTIFICATION 


By John C. Houbolt 


SUMMARY 


Treatment is given of system response evaluation, especially 
in application to subcritical flight and wind tunnel flutter 
testing of aircraft. An evaluation is made of various existing 
techniques, in conjunction with a companion survey report. 
Theoretical and analog experiments are made to study the identi- 
fication of system response characteristics. Various input 
excitations are considered. New techniques for analyzing 
response are explored, particularly in reference to the prevalent 
practical case where unwanted input noise is present, such as due 
to gusts or wind tunnel turbulence. Further developments are 
also made of system parameter identification techniques. 

Theory on the subject is extended, and many aspects of 
identifying system response characteristics are given in hand- 
book summary fashion. 


INTRODUCTION 


An important and vital phase of the aeroelastic study of 
aircraft is the substantiation of flutter by means of sub- 
critical flight flutter or wind tunnel tests. Because of the 
commonality of the problem to aircraft designs, a major confer- 
ence was held on the subject in May 1958 in Washington, D.C. 
Until recently, little had been done toward maintaining a 
summary of the various techniques used, or how they compare; 
essentially the companies have independently pursued and 
developed their own individual schemes. 

Communications, smaller group meetings, and conference 
papers indicate that there is much mutual interest in the 
subject, that a number of different methods are being used, and 
that a survey and critique of these methods would be valuable. 
Design considerations of the space shuttle system emphasize the 
need for and timeliness of such a survey. 

* 

In the aeroelastic analysis of the shuttle configuration 
many questions naturally arise. What subcritical flight flutter 
techniques are presently being used, and has the ever increasing 
modern computer developments led to improvements in the 
techniques? A significant question is whether the techniques 
used on aircraft are suitable for studying the flutter problem 
of the shuttle. It is realized that the flight parameters for 
the shuttle will be in a constantly changing state. Thus, 
steady-state-type flutter testing techniques may not be appli- 
cable, and those techniques based on transient excitation may be 
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the only type suitable. If present techniques do not appear 
suitable In application to the space shuttle, then research must 
be undertaken to develop flight flutter prediction methods which 
will be applicable. 

These needs led to the study effort that is covered in this 
report, sponsored by Langley Research Center of NASA. Part of 
the effort was directed toward making a survey, reference 1. 

This reference should be regarded as a companion to this report. 
Treatment herein deals with the evaluation of various sub- 
critical flutter testing techniques, with the set up and conduct 
of numerical and analog experiments of various schemes, with the 
development of improved procedures, especially for the case 
where input noise is present, such as due to turbulence or 
buffeting, and with the development of system identification 
techniques. The report is also intended to be, in part, a hand- 
book, since many notions used in system response evaluation are 
summarized . 

It is of interest to note that the survey and work of this 
report brought out the fact that not only is there much Interest 
in flight flutter testing in the United States, but a very deep 
rooted interest in England and other European countries as well, 
and that, in fact, several other survey-type papers on the sub- 
ject have recently been written, references 2, 3j and 4. 


SYMBOLS 


a a constant; elastic axis position from leading edge; 

lift curve slope 

A(co) real part of frequency response function 

b a constant; wing span 

B(a>) imaginary part of frequency response function 

c wing chord 

C(o)) amplitude of frequency response function 

e exponential function, position of c-g relative to the 

elastic axis, positive aft 

e „ distance of force application from elastic axis, 

positive forward 

e distance of accelerometer from elastic axis, positive 

forward 


f 


frequency in cps 
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F(t) 


Input force 


Fy(“) 


h 

h s 

HH 

k 


m 


L 

m 

n 

P 

q 

r 

r 


R y (t) 

S 

t 

V 

w 

y(t) 


Fourier transform of function y ; generally, the 
subscript denotes the function 

impulse response function 

sin as t 

response due to F r^- 

0) t 
o 

frequency response function, H = A + iB 
reduced frequency, k = ^ 
mass radius of gyration 
lift 

mass; an integer 
an integer 

operator ^ ; a root 
dynamic pressure, q = pV 2 
nondimensional elastic axis position, r = 

e 

nondimensional location of accelerometer, r = — 

* o c 

nondimensional location of input force, 
e f 

Y* = — - 

f c 

correlation function of y 

pvt, 

nondimensional time, a = =- L ^- 

* c 

time 

velocity 

nondimensional deflection, w = ^ 

c 

general displacement response function 

response due to noise 

deflection at accelerometer location 


0 damping coefficient 

P cr critical damping coefficient 

<5 Dirac delta function 


o|a> 


4 


€ 

X 

4 

V 


P 

a 


T 

0 

0y (a>) 

0) 


CO, 


CO 


time interval; also used to denote error 
characteristic root 
mass parameter 

2 v 

measure of velocity, v ~ — 
air density 

py 

reduced velocity, a - x___ 

r 

time 

angular displacement 

power spectrum of the function y 

circular frequency 

undamped frequency; a cut-off frequency 
reference frequency 


GENERAL THEORETICAL RELATIONS 


This section presents a listing of the principal general 
relations that are appropriate in the treatment of the response 
of linear systems to various forcing functions, and especially 
in application to suboritical flutter testing. For the most 
part, little associated discussion Is given. Some derivation 
is given where it is felt appropriate, especially where the 
relationships are not generally known or used. The equations 
are formulated in terms of basic concepts that are involved in 
describing the response characteristics of a structure, such as 
the impulse function, the frequency response function, corre- 
lation functions, and Fourier transform relations. Some of the 
equations presented may appear as new developments. 

General . - Let the general governing differential equation 
for response be given by 


V - V 


( 1 ) 


where Dq and D^, are differential operators, and y is the 
response to the forcing function F . For a simple damped mass 
oscillator equation (1) is 

my + 0y + Ky - F 


If the input force is a Dirac impulse function at t = 0 , 
equation (1) defines the impulse response function h for 
displacement as follows: 
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- D 2 6(0) 


( 2 ) 


For a unit sinusoidal input, F = e 


lout 


, and with 


y * He ia * 

equation (1) yields the frequency response function 

H(cu) = A(cjd) + iB(co) 
according to the equation 

(A 1 + iA 2 ) (A + iB) = N x + iN 2 


(3) 


W 


where A^_ , and A 2 , N 2 are the real and imaginary parts 
that are associated with the operators D^_ and D 2 . The A 
component of H is symmetrical with respect to the frequency 
a) , the B component is antisymmetrical . 

The h and H functions are related by the Fourier trans- 
form pair 


H 


-I 


he' 10 * dt 


(5) 




do) 


(6) 


Since h is zero for t < 0 , and because of the symmetry 
properties of A and B 3 it may be shown that the following 
relations also apply 


h 


•if 


A cos cut do3 


h 


a 

-? / 


B sin cut dou 


Useful limit properties of the H and h functions are 


H (0) = J h dt 

o 

00 

r \ 

h(0) = i J A dt 


(7) 


( 8 ) 
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• •• • 

The impulse response functions h and h , 
and for acceleration, follows as derivatives of 
defined by the relations 

X 

h = j' ico(A + iB ) e ia)t da) 


for velocity 
h , and are 



X 



-X 


X 

- f he~ iaA dt - A x + iB x 

— x 


h = - 




ou 

J oj 2 (A + iB)e icot d0) 


X 


H, 


‘ I 
00 

- 1 fc - 1 


TT lout , 
Hge dco 


-icut 


dt = A 2 + iB 2 


where 


A 1 - - ctB 
B 1 = oA 

2 

A^ = - a) A 
2 

B 2 = - o> B 


( 9 ) 


By the superposition theorem, the solution of equation (1) 
for any general forcing function F is given by 


X 

y = J F(x)h(t - t) dx 


( 10 ) 


The Fourier transform of this equation is 


p y ( tu ) = H(co)F F (aj) (11) 

where Fy and Fp denote the Fourier transforms, respectively, 
of y and F . Similar expressions in terms of h ( , h , , 

and H2 apply for the response variables y and y . If each 
side of this equation is multiplied by its complex conjugate, 
the result is 
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Vy = Hfl Vp 


which leads to the well known spectral equation 

( 12 ) 

relating the input spectrum 0p to the output spectrum 0y 
through the amplitude squared of the frequency response function. 


A specific response equation .- A specific form of 
equation (1) that is of prime concern in later sections of 
this report is 


tt t-tt .*» •• • ••• • 

y V + a^y iV + a~y + a g y + a x y + a x y = b^F + + b l F + b l F 

(13) 


With 


y 


.icut 

y. e 


(A + iB)e iajfc 


and 


t? ^ loot 
F = e 


this equation yields 

(A x + iA 2 )(A + iB) = N x + iN 2 


where 


A-l - oj 4 a 4 - o) 2 a 2 + a Q 
A 2 = or* - + oua^ 

N x = - A + b Q 
N 2 = - oPb^ + at» x 




Solution for A 


and B yields the results 
*1 A 1 


A 


N,A, + N p A 2 


b - -2-i 
B - 


2 “" 2 " 
A 1 + A 2 


N.^A n - N , A 0 


A 1 + A 2 
C 2 - A 2 + B 2 - 4^4 

A x + Ag 


1 ^ 
2" 


(14) 


(15) 
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Another solution of importance, but which has not grown 
Into popular use, lr. the solution for Impedance 
(Instead of for admittance A + IB). It is noted that 


1 

A +1'B 


A - iB 


From equations (15) and ( 9)3 it follows that 


A 

? 


B 


N 1 A 1 + N 2 A 2 
Nj + Ng 


N 2 A 1 ~ N 1 A 2 

+ n| 


- CD A „ 
L 2 

- od^B, 


( 16 ) 




C' 




1 


2 2 

A£ + Ag 

— 2 2 

N 1 + N 2 


i C 2 = A 2 + B 2 


These equations are of special significance in the treatment 
of flutter. 


A special type response function .- Let the input force be 
sin co r t 

taken as F = ^ — ; this Is a simple but very special input 

o 

force which has not received the recognition it deserves, and 
which gives a response that has many important and useful 
characteristics. For this force, equation (11) becomes 


F (cd) 


= H x 


-CD 


±JL- 

— /i\ 


CD, 


O 


CD, 


O 


where, symbolically, the notation indicates that Fy is simply 
H truncated by the "box car" function. The actual response, 
the Inverse transform of F , is herein denoted as h , and 

appears as y 


h 



sin uj t 
o 

CD T 
( ) 


h(t 


- t ) dx 


(17) 


Thus if the frequency response function of a system is known, 
and is sharply cut off to zero beyond a frequency cjdq , then 
the "impulse" function that is associated with the truncated 
frequency response function is simply the response of the 
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sin cu t 

system to an input force equal to — — j. 0 — : as or is made 

CO u U 

o 

larger, the more and more h_ approaches h . 

s 

Significant differential equations in terms of 
correlation functions .- Several equations are developed here 
which are not generally known but which should be of great help 
in system identification studies. 

Let the input force be h(- t) , which is the impulse 
response function folded around to fall along the negative t 
axis. By equation (10), the response would be 

00 

y •■= jf h(- x)h(t - t) dx 
-00 


which, in turn, may be written 

00 

y - j h(x)h(t + t) dx ( 18 ) 

-oo 

This equation is, however, to within some constant, the 
definition of the autocorrelation function of h(t) . Thus, 
the autocorrelation function of h is the response of the 
system to a force input of h(- t) , or 

D l R h = D 2 h <- ( x 9) 


A related equation is associated with the response of the 
system to pure white noise. Thus, if the input is white noise, 
equation ( 12 ) becomes 
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since 0 _ is flat. The inverse transform of this equation is 
r 


oo 

R ^ j h(x)h(t + x) dx 
y n -io 


which is the same as equation (18) . Thus, the correlation 
function of the response due to white noise is seen to be the 
same as the autocorrelation function of the impulse function h 


An equation involving the cross -correlation between any 
input and the associated response is also of significance. The 
Fourier transform of equation (1) is 


(A x + iAg)P y - (N 1 + iN 2 )P p 
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Multiply both sides by , the r omplex conjugate of Fp ; the 

result is 

(^1 + lAj.JPyPp - (Up + INgJPpPp 


which immediately leads to the spectral equation 


(A^ + ^2^Fy _ ^1 + ^2^F 


( 20 ) 


where 0p y is the cross-spectrum between the force and the 
response, and 0p is the spectrum of the input force. The 
inverse transform of this spectral equation is 


D R— 

1 Fy 


D 2 R F 


(21) 


Thus, if the autocorrelation function of a.n input F is applied 
to the system as an input force, the response is the cross- 
correlation function between F and the response y due to F . 


Some significant equalities and transforms .- Some 
equalities and transforms of importance are. listed here to close 
out this section. They are given without proof or derivation. 

Two basic Fourier transforms: 


y(t) 


. 1 


r 1 



-T 

T ♦ 


sin to t 
o 

v 




sin off 
off 



Convo lution : 
sin cd^t 

CO^T 


sin a> 2 (t - t) 
- t) 


sin oi. t 
d T =: X- 

Ol^t 

sin ciigt 
cu^t 

d. 

~ sin cii^t 
= i sin cu. t 

CL X 

= 0 


for 

<S) 2 

> 


for 


> 

“2 

for 

co 2 

> 


for 

“2 

= 


for 

cd 2 

< 

“l 


i 


sin 


<Y 


sin as (t - t) 

- t) dT 
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Properties of a specific h function: 


The impulse response function h for a simple damped mass 
oscillator is 

- — E_ oo t 

P,- o 


h =- 


1 „ ■'cr 


mco, 


sin co , t 
d 


( 22 ) 


where 


ou, = co. 
a o 




e 


cr 


The autocorrelation of this function, see equation (18), is 

Pc r 

cos co,t + — • — - - sin out 


"h = 6 

O B 
H cr 


— — CD t " 

Per o 


P' 


cr 


as reference 5 also shows. The derivative of R, is 

h 


Rt 


j P " fJU 0). 

o d B 
H cr 


e 


- § CD t 

®cr ° 


sin cd , t 
d 


d i 


(23) 


(2^) 


Interestingly, Rp is seen to be equal, within a constant, to 
h (equation 22)). It would be of interest to study what type 
systems have h functions that obey this property. 


CHARACTERISTICS OF VARIOUS INPUTS 


Spectral content .- Figure 1 indicates in summary fashion 
the spectral characteristics that are associated with various 
inputs that are of concern in subcritical flutter testing. 

Four distinctly different functions, the 6 function, 

sin a t 

white noise, a swept sine, and r-2— , are seen to lead to 

CD u 
O 

an ostensibly flat power spectra. Besides the 5 function, 
which is difficult to achieve in practice, the only function 
which leads to a truly flat power spectrum, and which extends 

sin oo t 

to zero frequency is the — function. The white noise 

(V 

spectrum is usually quite jagged. The spectrum for a swept 
sine has large lobes at the .Low and high frequency ends. The 
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sin cju t 

contrast between the swept sine and the — — -r — Is interesting; 

o 

the swept sine function has constant amplitude but varying f re- 
sin O) t 

quency; the — - - j. 0 function has a constant frequency but 

o sin oo t 

varying amplitude. The — - - x -- 0 function is a very attractive 

o^ 

function for use in system identification studies and has not 

been exploited sufficiently.' Because in the limit as oj -► » 
sin 0) t ° 

the — j.° function approaches the 6 function, it is 

o sin a) t 

suggested that — be termed the impulse sine function 

w o 

for ready identification purposes. 

Classification of the swept sine function .- The swept sine 
function has become rather popular for use in subcritical flutter 
testing, reference 1. The rate of sweep or total duration is 
one of the prime variables; with some tests the sweep rate is 
fast, in others the rate is quite slow.. For discussion and 
testing purposes, it appears desirable to make a classification 
of the rate or duration of sweep. The rate of change of fre- 
quency depends of course on the frequency range covered and the 
duration required to make the sweep. For the testing of most 
aircraft systems, however, it appears that classification can be 
based mainly on duration alone . The following classification 
is suggested: 

1) Fast sweep - one made with a duration of about 
5 seconds. 

2) Moderate sweep - duration of around 1 minute. 

3) Slow sweep - duration of around 5 minutes. 

Each of these sweeps has certain advantages, and certain 
deficiencies, depending on the application. Results later in 
the report will try to bring out some of the relative merits. 

Related differential equations .- It is perhaps of interest 
to note that differential equations associated with various 
swept sine wave laws, and with the impulse sine function, may 
be identified. This section shows, in the nature of an aside, 
the construction of these differential equations. 

Consider the differential equation solution to be of the 
following general form 

y = e* sin g(t ) (25) 

where f(t) and g(t) are any functions of interest. If the 
first and second derivative of equation (25) are formed, and if 
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these derivatives are combined in linear fashion with the function 
y , then the following, differential equation may be shown as a 
result 

y - ( 2 ? +4) .y +(- f + f 2 + g 2 + f I) y = 0 (26) 

v g v g 

Equation (25) is thus a solution of this differential equation. 

The differential equation for various swept sine laws and for 
the impulse sine function follow directly from these equations. 

Linear sweep law ; 

Consider that the solution is the swept sine wave often 
used as a forcing function in subcritical flutter testing, namely 


y = sin 6 = sin (a -f bt)t 


where instantaneous frequency is defined as 


CD - 


de 

dt 


or 


CD = a + 2b t 


(27) 


( 28 ) 


In terms of the beginning frequency cd 0 and the end frequency 
^ , after a sweep of T seconds, a and b are 


so that 


a = cd 
o 



1 

T 


CD --- 0) o + (0^ - CD o ) | 


Equation (25) reduces to equation (27) for 

f(t) ■== 0 

g(t) * (a + bt)t 

By equation (26), then, the differential equation yielding 
equation (25) as a solution is 

y - iHrW y + < a + 2M ’ 2 y = 0 


(29) 


Linear sweep down ; 
For this case 
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y = sin (a - bt)t (30) 

or f (t) = 0 , and g(t) = (a - bt)t . Equation (26) thus 
Indicates the associated differential equation to he 

y + & r ~ 2bt y + ( a - 2bt) 2 y = o (31) 

Exponential sweep : 

T at 

Let od = cd e 

o 

so that 

Q = g(t) = f cudt * ~ e at 

Through equation (26), with f(t) = 0 , the swept sine 

y - sin e at J (32) 

is thus found to be defined by the differential equation 

y - ay + o) 2 e 2at y = 0 (33) 


Linear period sweep : 
For this case 


CD ■= 


or 


and 


6 = g 


(t) -/ 


a - bt 
- 1 


y = sin 


cudt = — ^ log (a 


- ~ log (a - bt) 


- bt) 


(34) 


By equation (26) with f(t) = 0 , the differential equation is 
found to be 


• • 

y 


y + 



a - bt 


y = 0 


(35) 
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The impulse sine function ; 

With f(t) = - log out 

g(t) = tut 

equation (25) defines the impulse sine function 


y = 


sin out 

aSE 


( 36 ) 


By equation ( 2 . 6 ), the associated differential equation is found 
to he 


•• 2 • 2 

y+^y + aoy = o 


( 37 ) 


This equation is a special case of Bessel's differential 
equation, with the solution 

_ 1 

y = t 2 z 1 (- cut) 

2 


The function 


cos cut 
y out 


is also a solution. 


Other second-order time-varying systems : 

Equation (26), and the subcases given by equations (29), 

(31), (33), (35), and (37), are noted to be associated with 
linear systems with time-varying parameters. As a further 
aside, it may be noted that equation (26) may be useful in the 
study of various second-order systems having time -varying 
properties . 

A common approach in dealing with time -varying systems is 
to model the system and then to seek approximate solutions to 
the modeled system. Reference 6 is an excellent treatise along 
these lines. Consideration of equations (25) and (26) suggests 
an approach which is just the opposite. Thus, it is supposed 
that the solution is known; from the solution the differential 
equation is derived. This differential equation is then 
examined to see whether it represents the system being studied, or 
at least is a close approximation to the system. 

Consider, for example, f(t) = - pt and g(t) = cut ; 
equation (26) then indicates the well known damped oscillator 
equation 

y + 2py + (cu 2 + p 2 )y = 0 
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In I hi 3 case , a time invariant, second -order system is implied 
by the choices, tor f and g . By contrast, suppose 
f'(t) ■ at + bt : -’ , g(t) -at ; in this case, equation (26) 
indicates 

y - (2a + Ubt)y + 2b + (a + 2bt)‘ :1 + a^Jy = 0 

This equation is noted to apply to a system with a linear change 
in damping and a quadratic change in frequency. If, further, a 
and b are small relative to a , then the coefficient of y is 
roughly a constant. The equation would then represent a good 
approximation to a system with a linear change in damping. The 
nature of the homogeneous response behavior of such a system is 
in turn automatically given by equation (25). 


SECOND-ORDER SYSTEM UNDER CONSTANT FREQUENCY EXCITATION 


In studying the response characteristics of structures, one 
of the primary goals is to identify the frequencies and damping 
values of the various modes. Common or popular ways of identi- 
fying these quantities are summarized in this section in terms 
of a second-order system. 


Figure 2 depicts a popular type construction involving tf .. 
frequency response function, and specifically in the form of a. 
plot of ' B against A . This presentation is often referred to 
as the Kennedy-Paneu method, reference 7. Much discussion on 
this type construction is also given in reference 8. For a 
second-order system, the A's and B's are given by 


Displacement , 


F 

mar 

A + IB = 

1 - x + igx 


Velocity , 


lx 


mar 


A i + iB i 


: — 

1 - X, + Igx 


Acceleration, 


Ao 


+ .iB, 


I 


-x 


2 


T 

x 


F 

m 


+ igx 


d> , 

- ~ • “«• * 
o 


... p 


r 


where x 
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Characteristically, the plots for all three quantities resemble 
circles, and indeed the plot for velocity is a true circle. 

The resonant frequency is identified at the position on the 
"circle" where there is greatest arc length swept for equal 
frequency increments. Damping is found in two ways: the diameter 
of the circle is 1/g (assuming the response at zero frequency 
has been normalized to unity), or by the equation 


g = 


2Aa> 

oo tan 
o 


(38) 


where cu Q is the resonant frequency. 

The plot in figure 2d illustrates the results obtained if 
there is a mixture of viscous damping and structural "g" type 
damping; specifically, the frequency response function is 
given by p 

— 5 

mco 

A + iB * ^ 

1 - x + i(gx + g a ) 

All results shown in figure 2 are for g = 2 = .1 and 

g_ = .1 , and are presented on the assumption ££at the factors 
s 

-rs • BET » and I are 

mao o 

o 


Other means for evaluating frequency and damping are shown 
in figure 3 (again for a second-order system). The top sketch 
represents C 2 , the square of the amplitude of the frequency 
response function. The resonant frequency is associated with 
the peak of the curve . Damping may be found as shown, either 

from the peak value (if — ^ is known) or from the width at 

m*o 

the half -height position. As shown in the bottom of the figure, 
the h function, or the response that ensues after suddenly 
cutting off a resonant excitation, is still another way to 
estimate frequency and damping. The frequency is evaluated 
from the period T j damping is estimated from the decay of 
the peaks. The curves presented are based on viscous damping 
and provide a quicfc way for estimating damping from successive 
peak values. 


Figure 4 illustrates the impedance method which is another 
good way for estimating frequency and damping, although the 
schemes haven’t been pursued greatly. The plots represent 
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Damping and frequency are found as shown (again presented on 
the basis that mo^/F is unity) . 

Figure 5 is provided as a convenient reference figure to 
indicate the basic characteristics of the impulse response and 
frequency response functions for displacement, velocity, and 
acceleration for a second-order system. 


RESULTS FOR SIMPLE SYSTEMS WITH TIME -VARYING INPUTS 


The nature of the results that are obtained through use 
of swept sine and impulsive sine excitations are brought out 
in this section. Most of the results given have been obtained 
through study of a simple damped mass oscillator system, with 

an undamped frequency f = 9*95 cps , and a — = .05 • 


Figure 6 shows the results for a fast swept sine run, 
sweeping up from 4.8 cps to 24 cps in 4 seconds. The top 
figure is the input force, the second the response. The 
figure at lower left represents the autocorrelation of the 
input force; the function is seen to be composed of two 

s y type functions, thus implying a flat -type spectrum 

between an upper and a lower frequency. The figure on bottom 
right represents the autocorrelation function Ry of the 
response y . If the input force has a truly white spectrum, 
then R y can be shown to be the same as the autocorrelation 
function of the h function, see equation ( 23 ) and reference 5. 
Further, it is also found that, at least for a single degree of 
freedom system, frequency and damping indicated by the R^ 
function are the same as for the h function. It is of 
interest to note that processing in this case involves use of 
the response only; the establishment of Ry provides a ready 
means for estimating system frequency and damping. 


Figure 7 shows results for h as obtained from the y 
function of figure 6 by a randomdec-type technique (ref. 9). 

In contrast to the randomdec process described in reference 9> 
the process advocated here is developed in terms of functions 
that result from zero-crossings considerations; figure 7 (&) 
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Illustrates the zero-? rc.csingr, technique for constructing the 
randomdec signature . The thesis is that the randomdec signal 
so obtained is the impulse response function h . The result 
shown in figure 7(b) represents the sum of only 20 functions, 
starting at the point marked a in figure 6. As seen, the 
beginning portion represents quite well the decaying sine wave 
characteristic of the h function for a second-order system. 

The results at larger time values cannot be considered reliable 
because the randomdec summation involved only 20 terms. Here is 
a case where a moderate sweep rate would be better for randomdec 
purposes. If a moderate sweep had been U3ed, then the summation 
could have involved many more terms, with the consequence that 
h function derived would also be accurate at larger time values. 

Figure 8 presents results obtained by sweeping down from 
24 cps to 4.8 cps in H seconds. The autocorrelation function of 
the response y .is found to be virtually the same as for the 
sweep-up run of figure 6. 

Figure 9 presents the sweep-up results for the system with 
zero damping, and is given to show the caution that must be used 
in interpreting the autocorrelation function. The function 
shown :in figure 9 indicates that the system appears to have some 
damping. The response y , however, shows a persistence in 
oscillations, or ringing, after the resonance frequency is 
passed. This is a tip -off that the system has little or no 
damping. Here again, if a moderate sine sweep had been used, 
then many more cycler, oi persistence would be .indicated, which 
in turn would lead to a nondecaying autocorrelation function. 

Figure io applies to a fourth-order system wherein two 
modes have frequencies close together; in this case f ^ --- 9.yy 
and fo 10.73 cps . The response y appears as in figures 6 
and 8, but the auto< urrelation function of y Is seen to 
exhibit a beat pattern, as might be expected (see ref. 10). 

The high frequency and beat frequency seen can be used to 
estimate the system frequencies. The estimation of damping 
fox each of the modes, however, :Ls not easy. Techniques for 
handling the damping evaluation for such situations are in 
need of development. 

Figure 11 shows tie response that results due to use of 
an impulsive sine input. Shown at the bottom of the figure is 
the randomdec signature for h that is obtained from the 
response y (as outlined in fig. 7). The h obtained is 
virtually a true replica of the exact h for the system. 

This figure shows that the use of the impulsive sine function 
can be a. powerful tool for use In evaluating system response 
characteristics . 

The following table summarizes the frequencies and damping 
values that are indicated by the various evaluation techniques 
for the second-order system that was studied. 


REPRODUCIBILITY OF THE 
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cu 


Exact 

9.95 cps 

.05 

By R , sweep up 

.9.92 

.051 

Randomdec of y , sweep up 

9.92 

.050 

By R , sweep down 

9.92 

.051 

From h , impulsive sine 

9.96 

.048 


THEORETICAL FLUTTER MODEL 


As an aid In the study of various subcritical flutter 
testing techniques, a theoretical flutter reference model was 
developed. This model was used to provide exact answers; the 
model was also set up on an analog machine so as to provide a 
means for simulating subcritical testing. This section describes 
the theoretical model used. 


Differential Equation Formulation of 
Nonsteady Aerodynamic Forces 

A novel approach is given here for approximating the air 
forces that develop on an airfoil having nonsteady motion. The 
development automatically accounts for lag in lift effects, but 
avoids having to give explicit consideration to the commonly 
used F and G functions that are due to Theodorsen and 
Garrick for an oscillating airfoil, 

A study of oscillating airfoil theory and results for two- 
dimensional incompressible flow indicates that the basic lift 
forces on the airfoil may be represented as shown in the 
following sketch 
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where a refers to the position of the elastic axis of the air- 
foil system. Besides Lx and ]_,£ , an additional force and a 
moment associated with the inertia of the air act on the air- 
foil; these inertia forces will be neglected in this treatment 
and will be assumed to be taken into account simply as additions 
to the mass and rotary inertia of the airfoil. Expressions for 
and may be written 

t . , 

L x =• | pcV'b - y + c ( j h a ( t ‘ T ) dr (39) 


r apc^Vb i 
l 2 “ E 8 * 


(40) 


where h a refers to the lift that develops at the quarter chord 
position due to an impulsive change in angle of attack at the 
3/4~chord position, and b is airfoil span. 


The growth in lift on an airfoil following a step-function 
change in angle of attack is often given in approximation by an 
equation of the form 


1 


0(f) = 1 - 


V 


-V 


(f : 


Sometimes two or more exponential terms are included, but fVu 
present purposes a single term is considered adequate. A good 
approximation for a wing of finite aspect ratio is, in fact, 
(see ref. li) 

1 - 0 ( t ) - 1 - ,6e"*3 s 


where s - 


‘Vt 

c 


The derivative of equation (4l) yields h , thus 


h. 


a 


-b t 

(1 - a 1 )6(0) + a^e 


(42) 


where 6(0) signifies a Dirac function at t - 0 . 

The substitution of this equation into equation (39) gives 


L x - JpcVb 


( i - a x/(4) » a x b x / y(^) e • d ? 

where Y - V0 - y + - 2^0 . The derivative of 


W) 


is 
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L 1 = | pcVb (l - a^Y - J > Y( T )e- bl(t ^dr + 

■ 00 

(44) 

If equation (43) is multiplied "by b x , and the result is added 

M. 1 


to equation 


integrals, is found 


the following simple result, void of any 


4 + b l L l * T> P<=Vb[A - + bpt 


(“5) 


The development of the equations for flutter (or any arbitrary 
morion) can now proceed on the basis of this equation, rather 
than through a formulation which involves the F and; G 
functions; nonsteady aerodynamic effects will automatically be 
taken Into account. Note, a similar development can be made 
with two exponential terms in equation (4l) ; equation (45) would 
then appear as a second-order differential equation for . 


Equations for Flutter Model 

Consider the aeroelastic system depicted by the following 
sketches 




V 


0 




The quantity a d.-noV'S the elastic axis position; e repre- 
sents the distance of the mass c-g from the elastic axis; 
is the distance of external force application from the elastic 
axis, and e u denotes the position of a response transducer, 
such as an accelerometer, relative to the elastic axis position . 

(jti the basis of equations (40) and (^5), 
motion of the system can be shown to be 


the equations for 
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M-P 2 + M-kJ 

-larp'' -ip 

r 

-1 
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-nrp" 

M- -f P c + r P + p. k 2 

c c ^ 

- r l 

2a 2 p 2 +2b 2 p 

p 

-2r 2 a 2 p^-(a 2 +2b 2 r 2 )p-b 2 

p + b 


w 

0 

u 


where 


and where 




S - wing area - cb 

The response quantities of interest herein are for dis- 
placement and acceleration at the pickup location , These 

quantities are defined as ° 
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o 

Z W + — 0 

• • *• o l '° 

Z - W + — i 0 

From equations (46) the solution for 
f orm 

V IV 

a c z + a,,z + a„z + a^.z 4- a,z + a z 
6 4 3 - l ° 


07) 

(46) 

z is found to take the 
• « • 

-•= b^F + b^F + b x F + b Q F 

( z+ 9) 


which is the basis for equation ( 13 ; presented earlier 
(equation ( 13 ) is the result of normalizing by the coefficient 
a5) , Functions (14) and frequency response solutions (15) anil 
(16) are thus applicable to equation (49). If the right-hand 
side is set equal to zero, and z is assumed to be , the 

following characteristic equation of the system is found 

5 4 3 o 

A — a r \- + Si.X + a, X' + a„A 4* a. A + a — 0 
3 a 3 d i o 

The roots of this equation, A - 3 4 I to , are a key to the 
system response characteristics. El A Is taken as ia> , the 
characteristic equation yields the real and imaginary parts 

4- a 

o 

•f coa^ 

as defined also hy equations (.14). 

Equations (4b) (or equation (49)) were evaluated for a 
number of specific cases to provide exact reference results.. 
The equations were also set up on an analog machine so as to 
provide a means for simulating cube /it leal flutter testing. 
Some of the exact results are presented in the following 

i ater . 

by Flutter Model 
given herein are based on 


section; analog results will appear 

Theoretical Results indicated 

All of the theoretical results 
the following parameter choice 

q - 10 

" " L 

c" 


. a 2 

A^ a) a /( - co a 

c. o 

A . -• or 'a, - to* a,. 

6 3 
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r l - - 1 


' 2 




b, =• .3 


In all but one case, the uncoupled torsion frequency was 

taken as 10 cps, the uncoupled bending frequency coy was taken 
as 2 cps; in the lone case cu - cu The variable parameters 


were as follows 


a = 


JL ... 2V 

k CD c 
r r 


(cu^ is a reference freq.) 
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In all the result.’, to follow, the tel ere nee frequency cd. wa.. 
taken as 10 cps. 1 


System roots .- A common way to evaluate flutter speed and 
frequency is to solve lor the roots od of 

A x - 0 


A 2 = 0 

as a function of airspeed V ; the roots of Aq and A^ are 
referred to here as quasi-flutter roots, since they are 
fictitious values for all speeds except the flutter speed. A 
flutter condition is defined when the roots of Aq and Ay 
are equal. Figure 12 shows illustrative behavior patterns for 

the roots of A^ and A> due to variations in . Note tliat , 

in agreement with resylcs shown by other investigations, no 
flutter is found with the mass <■ „g. in a forward position, 

& 

~ - - .1 . It is also to be noted that if a sixth-order system 

instead of a fifth -order system were being studied, then the 
branch labeled a 0 would become: the. A^ branch. 

Figure 13 is for the lone case where ODy * uu. • this 
situation Leads to a very low flutter speed, as has often been 
observed. 

PAGC R 
^ Won QUALITY 





Illustrations of the behavior patterns of the true root:., 
as obtained from the characteristic equation, are shown in 
figure 14. These roots are of physical significance since 
they indicate the damping and frequency of the various modes 
that are present in any response excitation. Flutter occurs 
when one of the damping values (3) becomes zero (crosses from 
positive to negative damping) . Figures 14(c) through (e) are 
examples of a slow approach to flutter, since the damping de- 
grades to zero in a slow fashiop as airspeed is increased. For 
such cases a flutter speed prediction can usually be made by 
extrapolating the damping results. Figures 14(a) and (b) 
illustrate the behavior for an abrupt or explosive type of 
flutter. The damping may appear wel? behaved, but then with a 
very small speed increase can sudden Ly degrade to a positive 
value. These cases are very difficult, or impossible, to pre- 
dict in practice, and are the cases that cause grave concern in 
flight flutter and wind tunnel testing. A comparison of 
figures 14(a) and (b) is interesting. In one case, the fre- 
quencies of two modes cross, while the corresponding damping 
values diverge; in the other case, the damping values cross, 
while the frequencies approach one another but then diverge. 

For the other cases shown in figure 14, the frequencies tend 
to come together, but no crossing is noted. 

Figure 15 shows a comparison of the true and quasi-roots 
for frequencies. Figure 15(a) shows that for a mild approach 
to flutter (figure 14(c), the true roots and quasi-roots are 
markedly different (a true and quasi-root are of course the 
same at flutter). Figures 15(b) and (c), which apply to an 
explosive-type flutter, as seen in figures 14(a) and (b), show 
that one branch of the quasi-roots is close to the true roots. 
This fact, and other characteristics that are seen in 
figures 15 ( 1 ) and (c), as contrasted to figure 15(a), suggest 
that perhaps there may be a "tip-off" at subcritical speeds as 
to whether or not the flutter may be explosive. At least three 
distinct patterns may be noted. 

1) A quasi-root branch falls close to the true roots. 

2) The top portions of the quasi-root branches tend to 
remain parallel over a large speed range . 

3) The quasi-root branches tend to coalesce near one 
another (the coalescence points aren’t separated 
greatly along the x-axis direction) . 

Whether all these three characteristics must appear simultane- 
ously, or whether any one is sufficient to indicate the likeli- 
hood of an explosive flutter, is not known. Further study to 
investigate these tentative observations is considered desirable. 

Coefficients of the governing differential equation . - There 
is another promising possibility for extrapolating results for- 
ward to estimate the flutter speed, which should apply whether 



thu approach to flutter .in mild or of an explosive type. The 
scheme is based on the notion of identifying system parameters 
in a more complete way than identifying damping and frequency 
values only. Specifically, one concept is to identify the co- 
efficients a n in the governing differential equation of motion, 
equation (49). The problem may be stated as follows. Suppose 
the response z due to a given excitation force F is 
established; is it then possible to use z and F to estimate 
the coefficients a n , thus establishing the differential 
equation. If the coefficients are known, then the complete 
response characteristics can be determined. Involved also Is 
the notion that perhaps there, Is a more orderly variation of the 
coefficients a n with air speed than found for the damping or 
frequency values. Thus, the concept advanced is that of 
identifying the coefficients a n at several subcrltica.1 speeds 
assuming a certain order model applies, and then extrapolating 
these coefficients to higher speeds, and in turn to use the 
extrapolated values to predict the flutter speed. As a way of 
gaining some insight to this concept, a specific evaluation of 
the coefficients in equation (49) was made through use of 
equations (46) . Airspeed and elastic axis location were left 
as variables, the other parameters were given the specific 
values indicated previously in this section. - The results found 
are as follows: 

Sir-oP = 10 - 100r 2 
D t 


a ii h p 

* 5.72 - 30r - 8.6r 

= 10.4 + (.66 - 9 * 3 r)o 2 

= 4.3808 - (.15 + 3 r)a 2 
a l co r = + • 590 kc 2 


( 50 ) 


a 


o 


a 


.1 - .012a 2 


2V 

where a = . The coefficients are noted to be invariant 

r c 

or have a simple quadratic variation with respect to the air- 
speed. Figure 16 shows the variation of the coefficients, 

normalized to make unity, for r = ~ = .1 and ^ = .3 . 

Since the variation with speed is orderly, and since the 
theoretical model indicates the type of variation that each 
coefficient should exhibit (flat or parabolic), reliable 
extrapolation to higher airspeeds ought to be possible . The 
question is: "How well can the coefficients be evaluated from 


measured response data?" The circled points on figure 16(a) 
represent estimations from response data, and will be dis- 
cussed further in a subsequent section. 

Results of the form presented in figures 16(a) and 16(b) 
may have a special significance in predicting what type of 
flutter may be encountered, that is, whether flutter will be 
of the mild type or of the explosive type. Part (a) applies 
to a mild type flutter, see curve of figure 14(d). The 

variation of the coefficient with airspeed is seen to be 
slight. By contrast, part (b), which applies to an explosive 
flutter case (see curve in figure 14(a)), shows a much 

greater change of the coefficients with airspeed, especially 
the ao and ag coefficients. Marked changes in the 
coefficients therefore appear to be a clue or a "tip-off" 
that explosive-type flutter can be expected. 


Frequency response results .- Figure 17 presents repre- 
sentative results for frequency response as obtained from the 
theoretical model. Some results that were obtained from the 
analog computer are also indicated. The following table serves 


to show the 

parameters 

that apply 

to 

each figure: 

Fig. 

v = 2V 

£ 

ff 

e 

o 

c 

c 

c 

c 

17(a) 

60 

.1 

-.1 

.n 

17(b) 

100 

.1 

-.1 

. / Analog results 
also 

17(c) 

100 

.1 

.1 

.1 

17(d) 

100 

.1 

.1 

*3 

17(e) 

100 

• 3 

.1 

.1 

The general 

Intent is 

to show the 

variation in the frequency 


response function as brought about by changes in air speed, 
elastic axis position, and in the locations of the applied ex- 
citation force and the measuring transducer. 


Figures 17(a) and (b) also include the results that were 
obtained from the analog computer set-up of the system. The 
results indicate that the analog system duplicates the theo- 
retically exact equations of motion quite accurately. 


Figure 18 shows the nature of the plots that are obtained- 

A R 

when use is made of the impedance values and ~ , rather 

C 

than A and B . The contrast with figure 4, which applies to 
a simple second-order system, is to be noted. The loops in the 
impedance loci plots are associated with the higher frequency 
tildes of the higher order flutter system. Damping is measured 
oy the closeness of tne loops to the origin. Flutter occur a 
when a loop passes through the origin, thereby indicated zero 


damping . 
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Figure 19 ( a ) presents impedance results in the form of 

(see equations (16)) versus frequency. The dips or valleys in 
the curve are of special signficance, since the horizontal 
position indicates a mode resonant frequency, while the distance 
of the dip from the horizontal axis indicates the damping of the 
mode. The right side of the figure indicates the manner of 
using these valleys to extrapolate to the flutter speed. In 
this case, an extrapolated value precisely the same as the exact 
value is indicated. This figure is considered a significant type 
plot, since it gives a fairly complete picture of the development 
of a flutter mode, allows both damping and frequency to be t racket, 
readily, and leads to a fairly direct extrapolation to predict 
flutter. Figure 19(a) applies to a "mild" flutter case. The 
question naturally arises as to whether a similar type extrapo- 
lation procedure would be valid for an "explosive "-type flutter 
situation. In figure 19(b) results are given for an explosive 
case, specifically, the case covered by figure 14(a) . The marked 
curvature of the curve is perhaps a clue that the flutter may be 
of the explosive type. It is seen that the variation of the 

i curve is not as gradual as for the mild flutter case of 
C 2 

figure 19(a); the variation is not nearly as abrupt, however, as 
is noticed for the p 2 damping curve of figure 14(a) . Thus, 
this extrapolation procedure shows promise of applying to the 
explosive cases as well as the mild flutter cases. 


sin a) t 

A — ■ ,° VANE FORCE GENERATOR 
out 
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Oscillating vanes attached to the wing structure are used 
quite often as a means for generating an excitation force for 
subcritical flutter testing. The vanes are usually driven in a 
swept sine fashion to produce a swept sine force . By means of 
the simple form offered by equation (4-5), it is possible to de- 
rive an equation for vane motion which allows the generation of 
various prescribed excitation forces. 

Consider equations (39) and (4-0) in application to a vane 
executing rotary motion only about some axis (see sketch pre- 
ceding equation (39) ). The total lift on the vane is set equal 
to the desired vane force F , thus 

L = L x + L 2 = F 

From this equation the following relation may be derived 


f«l + b 1 L 1 + Lg + b-j^Lg = F + b-jF 

Thmiicrh mpsns nf pnna.ti nns f n.rui f 4 r^ . a.nH -hhp pniiflt.1 nn fnr 

O "X \ ~ / \ * - T. — — ~ 

Y following equation (43) (with the y motion suppressed), 
this equation may be written 





?L r 


*)* 


0 + 


b 1 0 


1 

27rqS 


(F + L X F) 


(^ 1 ) 


Solution of this equation for 0 for a stipulated F yields 
the vane rotational motion that is necessary to produce F . 
Equation (51) thus leads to a ready means for generating an 
impulsive sine excitation. In this case F is set equal to 
sin 03 t 

— • the associated solution describes the vane motion 

03 t 3 

o 

that is required to produce a vane excitation force of 
sin 03 t 

r-2— . it is noted that this development applies for the 

CO u 

o 

case of a fixed wing; application to the case of a flexible 
or movable wing should be satisfactory, however, as long as 
the force being applied to the wing is measured. 


DEDUCTION OF SYSTEM RESPONSE CHARACTERISTICS 


This section outlines the development of various tech- 
niques for identifying the basic response eharaeterist ic s of 
a subcritical flutter system. Attention is focused main 1 y ' r 
the frequency response function H , and the impulse respond 
function h . A primary objective of forced excitation 
testing is to derive such functions accurately so that 

a) damping and frequencies can be identified reliably, 

b) or that a much more detailed system identification 
may be made, such as the reliable evaluation of the 
coefficients of the governing equation of motion. 

One of the biggest problems of flutter testing is that of 
coping with an unknown "noise" input, such as due to atmospheric 
turbulence in flight flutter testing, or due to tunnel noise in 
wind tunnel testing. Since noise excitation represents such a 
serious obstacle in deducing accurate and reliable response 
characteristics, considerable attention was directed towards 
developing means for eliminating or masking noise effects. For 
this effort, use was made of an analog computer in conjunction 
with various input force generators and system response 
analyzers. This analog simulation proved invaluable, since many 
different schemes could be investigated rather quickly and 
repeatedly. 

Surprisingly, a number of different ways for coping with 
the noise problem were found, each having different merits. 

In the development of the procedures, the following inter- 
related questions were kept in mind (it is suggested that these 
points also be kept in mind as the results are read) . 
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1. What type of input is required? 

2. What record durations are required? 

3. Can reliable results be obtained from a single record? 

4. Should results be obtained in the form of a single 
record of relatively long duration, or should the 
analysis be based on numerous records of short 
duration? 

5. How much time is required to analyze the results? 

Point 5 is, of course, of prime concern in flight flutter and 
wind tunnel testing. Shut-down time, or ground time, between 
test points to await data analyses is not considered desirable 
in general. The desire is to make a test, analyze the results 
in a matter of seconds or minutes so that tests can proceed 
almost immediately to the next point. 

The following section describes the analog simulation 
system that was used. In the subsequent sections results 
obtained from the analog simulation study are described. The 
ideal situation of no noise in the input is described first. 

The following sections then describe 11 different schemes that 
were developed to eliminate, or at least minimize, the noise 
problem; 4 deal with discrete-frequency testing, 5 with time- 
varying excitation, and 2 deal with the use of response infor- 
mation alone due to noise (ironically, noise response can be 
used to establish noise-free response characteristics). 


Analog Set-Up and Associated Excitation and 
Measuring Equipment 

This section gives a brief description of the analog set- 
up that was used to simulate subcritical flutter testing of 
an airplane; whether the tests be in flight or in a wind tunnel. 

Figure 20(a) is a block diagram indicating the various 
pieces of equipment used. Four types of excitation input were 
used, a sine wave generator, a swept sine generator, a tape re- 
corder which could supply any specified input such as the 

~a>t ' ^" function, and a noise generator. In general, the noise 

input was treated as an unknown (nonmeasurable) quantity. The 
analog system was a representation of equations (46). Low-pass 
filters were used on both the input and output signals to ensure 
that the spectral content of the signals did not go beyond a 
certain frequency (to avoid aliasing). The Fourier analyzer was 
used as a ready means for processing the signals. The scope 
display and strip charts allowed for a "quick look" data analysis. 
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The printer and plotter allowed for the recording of data. 
Figure 20 (b) shows the analog schematic that was derived from 
equations ( 46 ), and which was vised for wiring the analog 
computer . 

Figure 21 is a picture of the simulated flutter testing 
system, set up at the Langley Research Center of NASA. While 
equipment from a number of different sources would be suitable 
for use, the following listing indicates the specific hardware 
used. 


EAI PACE TR-48 Analog/hybrid computer 

Hewlett-Packard 5^51 Fourier Analyzer 

5466A Analog to digital converter 
54-75A Control unit 
2100A Computer 
H01-3722A Noise generator 
5460A Display unit 
H51-18A Oscilloscope 
ASR 35 Teletype Corp. printer 
7046A x-y recorder 
5323A Automatic counter 
3403C True rms voltmeter 
Sangamo Sabre III tape recorder 

Tektronix R564B Storage oscilloscope 

Spectral Dynamics SD112 Voltmeter log converter 

SD127 MZ/TFA control 
SD104A-5 Sweep oscillator 
SD109B Co/Quad analyzer 
SD105B Amplitude servo/monitor 
SD122 Tracking filter 
Rockland 1022F Dual Hi/Lo filter 

Datagraph 5-510 CEC Strip chart recorder 

1-5H CEO D.C. preamp. 


The following sections deal with results obtained from the 
analog simulation just described. Most of the results given 
are for the following choice in parameters 


v 



100 



Cases which depart from these values are so indicated. Because 
of the limited storage capacity of the Fourier analyzer used, 
most of the runs involving the use of the Fourier analyzer were 
made covering a duration of 5 seconds only. Swept sine studies 
were therefore restricted to fast sweeps only. 
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THE 


Noise -Free: Inputs 

Four techniques are described here for establishing H or 
h when the input is free of noise. 

Dwell .- Figure 22 illustrates the frequency response re- 
sults that were obtained by a frequency dwell technique. In 
this case, the excitation force is set at a certain level and 
at specified frequencies. The response and input force are 
analyzed jointly by means of a Co-Quad analyzer to yield A , 
the component of the response in phase with the sinusoidal 
excitation force, and B , the component 90 degrees out of 
phase. This is a good technique if the time of dwell at each 
frequency is not a limiting factor. 

Swept sine input with Co-Quad analyzer . - Figure 23 is 
typical of the results that are obtained “by use of a swept 
since force input and the use of a Co-Quad analyzer to es- 
tablish A and B as in the dwell case. The main trouble 
with this method is that results depend on the sweep rate, 
and whether the sweep .is up or down. Generally, a sweep up 
tends to distort peaks towards the right of the correct value, 
while a sweep down distorts peaks to the left. Damping indi- 
cated is higher than actual. To avoid distortion, a very slow 
sweep rate must be used, and in such cases the dwell technique 
might just as well be used. 

Fourier transform approach using a swept sine input . - The 
basis for this approach is equation (11). A swept sine input, 
which covers the frequency range of interest, is used for 
excitation purposes. The frequency response function is then 
evaluated from the Fourier transforms of F and y according 
to the equation - 

II -- / (52) 

r F 

Results generally do not depend on the rate of sweep. 

Figure 24 illustat.es some analog results that were deduced 
by this procedure. The top sketch applies to H 2 , and h , as 
might be obtained through acceleration measurements. The second 
sketch from the top applies to [1 and h , as would be obtained 
from strain gage measurements. The sketch near the r center 
illustrates a novel way to evaluate the function ru 2 . If the 

functions H Q and H are perfectly formed, then the ratio 

h 2 2 

should evaluate to of , see equations ( 9 ). The sketch 

shows the rat^o as obtained from the analog results; the near- 
ness to an of variation is not bad, considering that no 
attempt was made to establish H 2 and H as accurately as 
possible. Both Ho and R are "'of Interest in practice; the 
lower frequency modes tena to be emphasized by the H 
function, while the higher frequency modes are emphasized by 

H 2 . The bottom sketch or, the figure applies to and h . 


REPRODUCIBILITY of 
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Impulsive slue excitation.- Figure 25 shows the nature of 
the response that was obtained from an impulsive sine excitation. 
Figure 26 shows the h s function that was obtained for a 
slightly impure impulsive sine force input; also shown is the 
associated frequency response function. The rather sharp cut- 
off of the function due to the use of an impulsive sine input 
is noted. 


Discrete Frequency Testing with 
Noise in the Input 

Frequency dwell.- Figure 2 7 shows results obtained using 
the Co -Quad analyze r . Noise causes the A (Co) and Bo 
(Quad) evaluations to be time varying. Results from tests at 
discrete frequencies thus lead to vertical marks, or inkblots, 
as shown. Time of dwell at each frequency was from 30 to 60 
seconds. Statistical 1 y , it seems reasonable that the centroid 
of each ink blot should correspond to the noise -free situation. 
Results for no noise are indicated by horizontal tics and, 
indeed, these tics appear associated with the centroid location. 

The curves on the figure are the result of sweep up and 
sweep down tests , using a Co-Quad analyzer to establish A 2 wju 
B p . It is seen that , with input noise present, such tests are 
of little value . 


Figure 28 shows results obtained by using different 
averaging times in the Co-Quad analyzer. The curves represent 
transfer loci obtained by sweeping up arid sweeping down wi. bh no 
noise present. Th*- "random walk" patterns shown below are cue 
to a noisy input md. represents results that would appear at 
the specific lo< rations A arid I', of the trammer locus if 
frequency dwells w '■oncucted at the frequencies indicated . 

It is seen that in-, r cased averaging t ime greatly reduces the 
size of the random walk pattern obtained; these results thus 
show the effect of using longer record lengths to "average out' 
noise effects. The consequence, of course, is a much increased 
testing time. Errors due to sweep rate and averaging time of 
the Co -Quad analyze r are treated in reference T.i. 


Correlation of i nput and outp ut.- The top sketch of figure 
29 denotes a discrete’ frequency ICnput excitation. The second 
sketch is the response recorded; it is composed of the response 
to F as well, as the response due to the unknown noise ex- 
citation. If the noise and applied force are uncorrelated, 
which they should be in this case, then the technique of formu- 
lating the cross -correlation function between the input force 
F and the response y should lead to noise elimination. The 
bottom sketch is the cross-correlation function found. It is 
seen to be quid: clear and regular* . A comparison of ftp with 

F in turn 0.1:1 own for the accurate evaluation of A and B for 
the 1. roquciu ,v -'oil idered . Those results are for only 5 seconds 
of record. Lon go 1 time dur.-xt Ions would lead to even better 


results 1 'of 
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Peak shifting .- In figure 30, F denotes a discrete fre- 
quency input force ; y represents the measured response. If 
the sinusoidal curves were not present on the figure, y would 
seem to be a response to noise only; it contains, however, a 
definite sinusoidal component. The concept of using peak 
shifting or peak enhancement can be stated with reference to 
F . Consider first the trace as given, next consider a like F 
function and shift it so as to align peak b with a , consider 
another like F ‘ function and shift to make peak c align with 
a , and so on; then add all the results. The result is the 
trace labeled 2 F . Note, shifting on the Fourier analyzer 
system is such to cause the information that is shifted off the 
left side to spill around in belt fashion and appear on the 
right side. Summations in the overlap region are thus not valid; 
because of this overlap problem, the ends of the 2 trace have 
been cleared. Next, do the same operations with the y function, 
using precisely the same shifts as for the F function. The 
result is 2 y . Frequency response information can then be 
evaluated from the 2 F and 2 y functions. 

The concept in this technique is that the shifting and 
adding operations causes the meaningful or intelligent part of 
the record to be enhanced, amplified, or reinforced, while the 
noise level remains the same. Note, the addition of a number 
of uncorrelated noise records should give a result which Is 
similar to any one record. Thus, if an average value of the 
summation is formed, then the peak level should be the same as 
the original peak level, but the noise content should be de- 
creasing as ^ , where n is the number of samples involved. 

The establishment of a fairly clean sine wave, as indicated 
by the 2 y , from the rather noisy function y indicates the 
usefulness of this peak shifting technique. The result is for 
only 12 additions; a larger number would lead to an improved 
quality for the 2 y . 

Ensemble averaging .- This technique is based on the concept 
of adding together a number of independently generated records, 
with the additions being made so that the input records start 
in the same way (phase maintained) . As in the peak shifting 
technique, the ensemble sum of the output should show a de- 
creasing noise content as the number in the ensemble is in- 
creased. Figure 31 illustrates the procedure. The top sketch 
denotes a single input record, the second sketch the measured 
response. The records at the bottom shows the ensemble sum for 
F and y as obtained from an ensemble size of 30. Note how 
the noise in the output has been reduced. Frequency response 
information follows from the 2 F and 2 y records. 

Time-Varying Excitation with Noise in the Input 

The use of a swept sine input and the Fourier transform 
relation given by equation (52) leads to results of the type 
shown in figure 32 when an unknown noise is present in the 
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input. This figure should be compared with the no-noise analog 
result, figures 22 and 26 , and with the corresponding exact re- 
sult, figure 17(c). This figure vividly illustrates the problem 
brought about by input noise. With results of the type shown, 
it is virtually hopeless to deduce meaningful response character- 
istics. Techniques for obviating the noise problem are there- 
fore of great interest. 

Clearing or weighting of the h function .- It is to be 
noted that throughout this report, discussions of H or h are 
essentially synonymous, since, as equations ( 5 ) and ( 6 ) show, 
knowledge of one function automatically defines the other. The 
use of equation (52) leads first to H , but h then follows 
directly. Certain features in figure 32 are worth noting. The 
position labeled a on # h appears to be the point where the 
correct or noise-free K function would have decayed to near 
zero. Beyond this point the information shown is mostly due to 
noise. In turn, most of the jaggedness in the kg and Bp 
functions is due to the erratic behavior of h beyond a point 
such as a . The simple technique of clearing the h function 
beyond the point a is thus suggested as an easy means for 
vastly improving the noise problem, reference 13. A rectangular- 
truncation is implied, having unity out to a selected time, and 
zero thereafter. Figure 33 illustrates results of this type 
clearing process. The results on top is another example of the 
type of results shown in figure 32 , and are for the raw data. 
Clearing the h function beyond a point corresponding to a on 
figure 32, and then rederiving H leads to the results shown in 
the middle. A great improvement is noted; noise effects are still 
present, but at least some indications of frequencies and damping 
are present. The bottom sketch applies to the no-noise case and 
is included for comparative purposes. 

References 3 and 13 indicate the use of an exponential 
weighting function on the raw h function, and then re- 
evaluating H , as a means for minimizing noise. Figure 3 A 
shows results obtained by this approach. This technique appears 
to be quite effective in eliminating noise effects. The intro- 
duction of a weighting function of course causes distortions m 
the derived H function. Corrections that account for the 
weighting function must subsequently be made to the deduced 
damping values. Fortunately, these corrections are easy to 
make when an exponential function is used, since the correction 
is simply to subtract out the apparent damping that has been 
added by the weighting function. Other weighting functions are 
not recommended, since the corrections are not known or cannot 
be made . 

Cross-spectrum between F and y .- The theory for this 
technique is as follows.' The response due to an applied ex- 
citation and due to the noise environment may be written as 

= /(F+ F n )h(t - T)dT 


y = y p + y n 


(53) 
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where Fn is the unknown input noise and y n is the associated 
noise contamination in the response. The Fourier transform of 
y is 

F = F + F = (F + F )H 
J y F n 



If this equation is multiplied through by the complex conjugate 

F , then the following spectrum equation is indicated 
F 

0 Fy F + 0 Fy n = + 0 FF n ^ H 


Because there is no correlation between F and F n , however, 
both the cross spectra 0 Fy and 0 FFn should vanish. The 
equation then yields n ^ 

tt . ^f 
0 p 

which appears as a completely noise-free result. The technique 
is thus to form the cross-spectrum 0 Fy between the applied 
input and the measured response, and to divide by the input 
spectrum to obtain H . 

Typical results are shown in figure 35 for the case of a 
swept sine excitation. A substantial improvement is noted. 

The illustration is not a fair test of the approach, however, 
because of the very short record lengths that had to be used. 
The results shown in figure 35 represent only 2.5 seconds of 
data, because correlation was involved. In spite of this 
limitation, a marked improvement in the noise problem is noted. 
It is felt that record lengths of about 30 seconds (a moderate 
sweep) are needed for this correlation technique, and that if 
such lengths were involved, then almost perfect results would 
be obtained. The author considers this to be one of the best 
techniques available for eliminating noise effects. 



Peak shifting.- A peak shifting technique similar to that 
described under discrete frequency testing with noise in the 
input is also possible for a swept sine input excitation. In 
this case, shifts are based on the peaks of the swept sine in- 
put function; the shifts for the response are taken identical 
to the shifts of the input. The sum of the input shifts is 
treated then as a single input function, and the sum of the 
shifted output functions is treated as a corresponding single 
response function. Note, the concept that solutions for linear 
systems may be linearly added is involved. The summed results 
are treated by equation (52) as though they represented a single 
response run. Results obtained are shown in figure 36. Rela- 
tive to figure 32, a vast improvement is found. Here again 
though, as with the correlation example, the test of the tech- 
nique is not fair. The short record length available allowed 
only a small number of peaks to be summed. Thus, enhancement 
of the meaningful signal part was not sufficient to "average 
out" the noise. Record lengths of around 30 seconds (moderate 
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sweeps) should allow construction with many peak shifts and 
should make this a powerful technique. 

Ensemble averaging .- The previous techniques are based on 
the use of a single record only. By contrast, an ensemble 
averaging technique involves the use of many sweep runs. The 
idea is simply to evaluate H or h for each of the runs and 
then to add the results to form an ensemble average. The con- 
cept is that noise effects will "average out" to zero. Results 
for an ensemble of 20 sweep runs, each of 5 seconds duration, 
are shown in figure 37. The figure gives the results for both 

• • • 0X7* 

h and h , and for values of v = — = 60, 80, 100, and 110 . 

The value of v for flutter is 117 . This technique is seen 
to be very effective in leading to reduced noise effects; the 
penalty is the problem of having to make many repeated runs. 

If 20 to 30 runs, each of 5 seconds duration are needed, 
then a total test time of 100 to 150 seconds is implied for 
each test point (one speed), not counting reset times between 
runs.. A total time of 150 seconds approaches but is still 
smaller than the sweep duration runs of 4 minutes that are 
often used in flight tests. The question is raised: "Would a 
single run of only about 50 seconds, analyzed by the corre- 
lation or peak shifting technique, be better?" Unfortunately, 
this question can not be answered at the moment. 

Figure 37(a) also includes the result for | | > which 

also represents the spectrum of h . This function is' seen 
to be quite clean, and of all the functions shown, allows for 
the easiest evaluation of system damping and frequency. The 
frequency is indicated by the location of the peak, the damping 
by the width at half -peak height, see figure 3* 

Combined ensemble averaging and h weighting .- Some of 
the methods described here can of course 'be used in combination. 
Figure 36 shows the results obtained by an ensemble averaging 
of only 5 runs, with the subsequent use of the exponential 
weighting function technique. The resulting curves are quite 
smooth, but correction of the data must of course be kept in 
mind; the results shown should be compared with the exact re- 
sults shown in figure 17(c). 

Ensemble averaging using response to noise only . - The 

treatment following equation (19) indicates that the auto- 
correlation function of the response to white noise is the 
same as the autocorrelation function of the impulse function 
h . This fact suggests that a useful result might be derivable 
by working with noise response records only. Autocorrelation 
functions of the response to noise alone were established for a 
number of individual runs. The result of adding together 20 
such functions is shown in figure 39. The agreement with the 
result shown in figure 37 is remarkable. Damping and frequency 
appear to be readily identifiable. The Fourier transform of 



39 


the autocorrelation function is shown as the second function 
from the top. The smoothness of this function indicates that 
the results are essentially noise-free response results for the 
system. The Fourier transform of the right half of Ry is 
often of interest. This result is shown at the bottonf of the 
figure . 

Use of the randomdec technique .- The randomdec technique 
is another means for deriving system response characteristics 
from noise response information only. The technique described 
in figure 7 was applied to the noise response of the flutter 
system. Results are shown in figure 40; these results are to 
be compared with the h results shown in figure 24. 

The randomdec technique may be applied to acceleration 
noise response results, but the construction does not lead to 
a system physical function. The reasons is that the Dirac 
function usually found with h functions, see figure 5* are 
not accounted for properly. The randomdec signature found 
from acceleration should give, however, an indication of 
system frequency and damping. 


SYSTEM PARAMETER IDENTIFICATION 


Most of the techniques used for evaluating frequency and 
damping of the various modes of a multimode system are based 
on the behavior of a single degree of freedom system. When the 
modal frequencies are well separated, reasonably good estimates 
of frequencies and damping probably result, but even in such 
cases, the values deduced are really only "pseudo" estimates of 
the true values. When frequencies are close together, identi- 
fication becomes uncertain or impossible, or estimates may be 
in large error. The identification of the parameters of a 
system in greater detail is therefore desirable. The establish- 
ment of the coefficients of the governing differential equation 
of motion (see equation (IJ)) from response measurements is, 
for example, an extended step to better system identification. 
With these coefficients, all response characteristics of the 
system can be evaluated, whether frequencies are close together 
or not. The accuracy of the determination of the coefficients 
is of course a factor. 

This section describes three means for evaluating system 
parameters in greater detail from response measurements. 
References 14 through 17 represent noteworthy treatments of the 
subject . 

Collocation using the frequency response function .- Assume 
that the system under consideration is a 5th-order system, so 
that equation (13) applies. The frequency response equation 
for displacement response indicated by this equation is 


RLPRUDUChji.ai i or i iixj 
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t 


co 4 a^ - a) 2 & 2 + a Q + i(oP - oPa + cue., )|(A + iB) =-- 




- a + b .. -I- 1 


o 


i(-oPb.~ + oib^ 


(56) 


which when expanded leads to 


co^Aa^ - co 2 Aa 2 + Aa Q + opBa^ - ofia^ + oPbp - t> 0 - oPb (57) 


4 

2. 

to Aa^ - 

co Aa 

4 

2 

a) Ba^ - 

co Ba 


co^Ba^ - 6o 2 Ba 2 + Ba Q - op Aa^ + + oPb^ " = _a} A (58) 

The multiplication through of equation (56) by A - iB and 
division By C 2 = A 2 + B 2 leads to the following two 
alternative forms 

m\ - m 2 a 2 + a o + m 3 Sj b 3 + m 2 ^ i 2 - cb 5_ ^ b Q = 0 

(59) 


- oP a + cua, + op t - oP ^ b - co ^ b, + ^ 

3 1 c 2 3 c 2 2 c' 1 C d 0 


(60) 


The collocation solution proceeds by using these equations 
singly or jointly to solve for the a n and b n coefficients. 
Consider equation (57) for example; five a n and two b n coef- 
ficients appear in this equation. Measured values of A and 
B are substituted in this equation at seven different values of 
oo , leading to seven linear simultaneous equations with unknowns 
a n and b n . Solution is then made for these coefficients. 

Or, five frequency values may be used in equation (57) and four 
in equation (58), giving nine simultaneous equations in terms of 
the total of nine unknown coefficients. 

A test of the approach was made by using values of A and 
B as obtained from the exact solution, equation (49) . The 
coefficients evaluated were in good agreement with the original 
coefficient used to obtain the A and B values; this compari- 
son is shown in figure 16 , where the circled points refer to 
the coefficients as evaluated by the collocation procedure . 

This comparison indicates that the scheme works, at least in 
principle . 

Analog values of A and B were also used to check the 
procedure. Some of the results obtained were good, some were 
bad, depending on the frequency locations chosen. Figure 41 
indicates a few of the results obtained, and shows the quasi- 
roots that were obtained from the deduced coefficients in 
comparison to the exact quasi-roots. The results indicate 
that perhaps the best procedure to use is to evaluate the 


/ 
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coefficients several times for different frequency choices (the 
evaluation is very quick since only a few simultaneous linear 
equations are involved) and then to average the results. 

To check further on the collocation procedure, a sensi- 
tivity study was made to establish how sensitive the coef- 
ficients a n and b n were to assumed changes in A and B . 
Figure 42 indicates the results of this study. The study 
started with the exact values cf A and B . These values 
were given random variations through use of a random number 
generator. With the varied values, solution was made for the 
a n and b n coefficients. This, experiment was repeated 100 
times. The value in figure 42 represents the standard 

deviation of all the variations of the A and B values; the 
value, the standard deviation of all the variations found 

values. This figure shows that the 
e sensitive to the A and B values 
used (roughly, a magnification in errors of 200) . 

Least squares difference equation approach . - The 
difference equation equivalent of equation ( 13 ) may "be written 
as 

y n+ 5 + a 4 y n+4 + a 3 y n+3 + a 2 y n+2 + a l y n+l + a o y n = 


for the a n and b 0 
coefficients are quit 


b 4 P n+4 + b .3 P n+3 + b 2 P n+2 + b l F n+l (o'.) 

where the y’s and F's represent equally spaced values with 
interval e . The a n and b n coefficients used here are not 
the same as the coefficients in equation ( 13 ), but rather are 
some combination of these coefficients. Assume the y n and 
F n values are measured values as obtained from a test, and re- 
write equation (6l) in the form 

y n+5 +a 4 y n+4 + a 3 y n+3 + a 2 y n+2 + a l y n+l + a o y n 

b 4 P n+4 “ b 3 P n+3 " b 2 P n+2 * ^l^n+l “ fc n 


where e n represents a possible error because the y n and 
F n values are not exact. The coefficients a n and b n are 
now found using a least squares process involving the error 
e n . The problem statement appears as 




0 


(63) 
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thus leading to nine linear simultaneous equations in nine un- 
knowns. Solution yields the desired values of a 

n 

With the solution for the a Q values, the roots of the 
char* teristic equation of the difference equation may be found. 
Assume the right-hand side of equation (6l) is zero and let 

y = e^ = e^ ne ; the result is 

e; 4. o p 

p + a^p + a^p^ + a ? p + a^p + a Q = 0 (64) 

where p = e Xc . Solution of this equation yields the roots 


p = c 

n 


± id 


n 


Roots which approximate the roots of the characteristic equation 
of equation (26) may be now established as follows. Let 
X = |3 - ico represent these approximate roots; then 


e 





(65) 


Solution for 6^ and a> n yields 


d 

tan cu n e = (66) 

n 


^^n € 2 2 

e n = c + d^ 
n n 


(67) 


.. Application of this least squares procedure was made to the 
h > and R h functions, as given in figure 37. Since these 

functions represent homogeneous solution of equation ( 13 ), no 
Fn terms had to be considered; that is, equations ( 63 ) were not 
involved. Note, the initial values of the functions were not 
used so as to avoid the problem of their Dirac function type 
beginning. Results obtained for and p for the value of 

v = 110 are shown in the following table. Trie estimates for 
the first mode show a considerable variation; second mode values, 
however, are in good confirmation. It should be kept in mind 
that these limited results do not represent a fair test of this 
difference equation approach. First, the spacing e was rather 
large; a value of e, ~ .02 sec. was used, and thus one cycle of 
the higher frequency mode is represented by only six points. 
Second, the maximum number of points used was only 64. And, 
thirdly, there is very little response of the first mode present 

•i • 

in the h functions. 
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Differential equation solution by leant squares.- An 
inte renting least ’3;quar - e:. .(TTublon can be made in the time plane, 
but which makes use oi frequency plane information. Consider 
equation ( 13 ) in the following form 


V 


+ a 


4 y 


IV 



b.,F - b 0 F - b,F 

j r - x 


b G P = e(t-) 

(68) 


where e represents an error due to the use of measured y and 
F values. Assume that the response data have been processed to 
lead to y ~ h , see equation (17) . The function h corre- 

O • i - o 

s j n cu t 

sponds to the response due to F = — — , and falls out in a 

o 

straightforward natural way wnen processing is made on the basis 
of the Fourier transform relation, equation (5^)* Note, this 
analysis, although presented in terms of h and 
sin cd t s 

F _ o_ no t restricted to these functions alone. 

CD t 1 

o 

The coefficients a ri and b n are now evaluated through 
use of a least squares statement involving e , thus 


E 



dE 

TET 

n 




e £ (t )dt - Min. 


-^2 


• 44 dt = 0 


- , n 
dt 


'1 


f . J ,n c , 

€ 4iL at 0 


dt 


n 


J i 


( 69 ) 


(70) 


(71) 


Equations (70) and (71) give simultaneous linear equations in 
terms of the a n and b n coefficients, which then allows their 
evaluation to be made. 


Equations (70) and (71) may be converted to a more amenable 
form as follows. With limius extending from t -- - oo to t - <* , 
the integrals which appear in the simultaneous equations may be 
written in three general forms 
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Use is now made of two significant Fourier transform theorems, 
namely 

,n 

.1 . The Fourier transform of is (ito) F , where F\, 

dt y 

is tiie Fourier transform of y 


s on 

F 1 (f)F 2 (- f ) df - ^.7-, (t) y ( 4- t ) dt where o> = 2-rrf 


With these theorems , equations 17') he come 
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For the y - h_ function, let 
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Equations (73) thus become simply 
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Solution of these equations leads to the differential equation 
coefficients a n and b n . The symmetry of equations ( 75 ), 
and repeated appearance of the D n , E n and F n coefficients 
should "be noted. In all, only 14 such coefficients need 
evaluation (the four F n values are known). One of the 
appealing merits of this approach is that the D n ’s and E n 's • 
may be evaluated automatically by the Fourier analyzer that is 
being used to analyze ’the response. Solution for the a n 
values should therefore be quite easy and quick. 

A minor test of equations (75) was made using only approxi- 
mate values for D n , En and F n ; results for a n were, 
however, in reasonable agreement with the corresponding exact 
values. Further study of the soundness of this least squares 
technique is worthwhile . 


CONCLUDING REMARKS 


Main emphasis in this report has been on the development 
of Improved subcritical flutter testing techniques, whether in 
flight or in the wind tunnel, and particularly in reference to 
the situation where input noise is present. Discussion was 
given to a certain extent as the material was presented. Some 
additional observations are made here in the form of concluding 
•remarks . 

Frequency and damping evaluation: 

A number of different schemes were outlined for deducing 
frequency and damping from response measurements. Mode fre- 
quency is identified quite well by most all of the techniques. 
Damping determination, however, is more nebulous. The transfer 
locus scheme, which hopefully leads to identifiable circles, is 
popularly quoted but the author does not favor this approach 
too highly. The diameter of a circle obtained from acceleration 

IF 

measurements is, for example, — — (assuming of course that 



4 b 


behavior is that of a single degree of freedom system). Thus, 
explicit g evaluation from the diameter is precluded because 

the unknown — is also involved (m , in general, is some 
m 

generalized mass value). Actually, it would be better to use 
equation (38) to evaluate g (if the center of the apparent 
circle can be fixed), and then to use the diameter to evaluate 
F 

— . The methods preferred are those shown in figure 3 » and 

particularly the scheme involving use of the width of the 
spectrum peak at half power. This scheme is simple and direct, 
and subjective interpretation is a minimum. (Note, with g 
established by the width of the peak, the height in turn may be 

used to estimate — .) The scheme of deducing damping from a 

randomdec signature is also considered good and reliable. It 
is to be noted that all of the schemes are subject to a common 
problem; specifically, all the methods for deducing frequencies 
and damping are open to question for the situation where two 
(or more) frequencies of the system are close together. Means 
of being able to detect when frequencies are close together and, 
in turn, of deducing the frequencies and damping are in need of 
further development. 

.Deduction of H or h : 

As indicated, a number of methods may be used for obviating 
the input noise problem. The ensemble averaging technique is 
attractive but requires a substantial number of repetitive runs. 
The use of exponential weighting of the h function appears 
satisfactory, but schemes which do not lead to distortion of Uu 
data, which then require correction, are judged preferable. ‘Inc- 
cross -correlation approach (equation (21) and equation (55)), 
and the peak shifting technique seem the most appealing on the. 
overall. The technique of using the response only, as obtained 
from a "white" noise environment, and forming the autocorrelation 
function, with ensemble averaging, is quite intriguing and should 
be used where possible; in this case, a function representing the 
autocorrelation function of h is found. In general, randomdec. 
signatures may be interpreted as h only for displacement 
response. Randomdec signatures obtained froip velocity or 
acceleration response are not strictly the h or h functions 
because step functions or.jDirac functions at the origin are not 
reproduced (see h and h in figure 5). Damping and frequencies 
as established from the randomdec signatures for velocity or 
acceleration response should, however, be representative of 
actual system response characteristics. 

System identification: 

A variety of possible system identification schemes were 
developed, but attention was restricted herein to three 
approaches. The collocation scheme is quite simple but suffers 
from the fact that the frequency values chosen for response 
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matching is so arbitrary. The question of how to handle the 
two equations in combination (one associated with the real 
part of the solution, one with the imaginary part) is some- 
what of a mystery. The difference equation approach is con- 
sidered good and quite attractive. An analogous scheme, 
reported in references 2 and 13, appears to be highly regarded. 
The differential equation approach, which makes use of fre- 
quency plane information, needs more corroborative study, but 
is considered quite promising. 

A final word is given with respect to subcritical flutter 
testing techniques which differ from the type discussed herein. 
Tacit in the schemes mentioned in this report is the assumption 
that the air density is constant, and that tests proceed on an 
incremental increase in air speed basis. Reference 17 describes 
procedures for testing on an increase of air density basis, 
holding speed essentially constant. Figure 43, taken from this 
reference, shows the excellent success that was obtained in 
extrapolating to flutter through use of a density increase 
approach. This technique is still considered to be a good and 
useful approach and should be kept in mind in any further 
development studies of subcritical flutter testing. 

Application to the space shuttle : 

Flight of the space shuttle will represent a situation of 
a time-varying system, since dynamic pressure and Mach number 
in particular will change rapidly with time, see figure 7, of 
reference 1. The. question that naturally arises is whether 
the methods discussed in this report, which apply mainly to 
time invariant systems, can be used for subcritical flutter 
evaluation of the shuttle system,. With respect to the various 
methods, the following recommendations are made. During the 
transonic region of flight (where the dynamic pressure will be 
nominally about 75# of q max ), the random forcing input associ- 
ated with the transonic flow will probably be strong enough to 
excite a sizable random response. It follows then that two of 
the methods discussed herein might be useful in evaluating the 
shuttle response characteristic. One is to break the record 
up into 5-second segments (on the assumption that the system 
is nearly time -invariant during such an interval) and then to 
form the autocorrelation function of the response for each seg- 
ment. The second approach is to form the randomdec signature 
from the random response record. 

At maximum dynamic pressure, the Mach number will nominally 
be around 1.5. It is not known at this time whether the flow 
during this period will be rough enough to cause random 
excitation. If tests indicate that the flow should still be 
"rough, " then the same two means for analyzing the data should 
be used. 
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To ensure that system excitation can be obtained over all 
regions of flight, it is desirable to have aerodynamic vanes 
or inertial shakers installed. The excitation recommended is 
to go through a sine sweep, first up, and then down, in 
continuous succession. Sweep durations of 5 seconds are 
suggested. Each five seconds of response information could 
then be analyzed to deduce H and h through means of the 
method involving the Fourier transform of the output response 
to the Fourier transform of the input force. In this case, 
the procedure involving the peak shifting technique would 
appear ideally suited. Note, the time-varying aspect more or 
less precludes any ensemble averaging approach. The auto- 
correlation of the response due to the sweeps should also be 
obtained as an additional means for evaluation of the response 
characteristics. 

In summary, three approaches appear useful for subcritical 
flutter test evaluation of the space shuttle system, namely 

1) The autocorrelation of system response, whether the 
response is due to a natural random excitation or 
due to a controlled force excitation. 

2) The randomdec signature approach. 

3) The use of the peak shifting approach. 
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Fig. 1.- Spectra for various input forces 
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Fig. 3.- Frequency and damping as obtained from the frequency 

response function or from the decaying free oscillation . 
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b) Randomdec signal for h obtained from y in figure 6 


Fig. 7.- Establishment of h by randomdect, technique 
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Fig. 12.- (Concluded) 
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ig. 14.- (Cont.) 
















Fig. Id.- Exact a n values compared with values deduced from frequency response function 
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Fig. 18 Illustrative impedance loci plots. 
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(b) ^ = .3 ; explosive flutter 


Fig. 19.- (Concluded) 
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Fig. 28.- Dwell results obtained with noisy input for various 
averaging times in Co-Quad analyzer. 























Fig. 34.- Improved frequency response function by exponential 
weighting of h function. 
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Pig. 37.- (Cont.) 
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of 20 runs 


Fig. 37.- (Cont.) 
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Fig. 39.- Ensemble averaging of the correlation functions for 
acceleration response due to noise input only, 
v - 100 . 
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Fig. 41.- Quasi-roots obtained from experimentally deduced differential equation 
coefficients. 
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Pig. 42.- Sensitivity study of deduced differential equation 
coefficients. 






